NBIS Report

01-Aug-2022

NBIS ID: SMS_6198
Report Version: 1.0
Request by: Mona N. Högberg (mona.n.hogberg@slu.se)
Principal Investigator: Mona N. Högberg (mona.n.hogberg@slu.se)
Organisation: SLU
NBIS Staff: Juliana Assis (juliana.assis@nbis.se)

1 Setup

## LIBRARIES
library("dada2")
library("devtools")
library("dplyr")
library("ggplot2")
library("microbiome")
library("phangorn") 
library("phyloseq") 
library("Rcpp")
library("reshape2")
library("tidyr")
library("vegan")
library("ShortRead")
library("Biostrings")
library("DECIPHER")
library("SensusR")
library("gplots")
library("gridExtra")
library("grid")
library("ggpubr")
library("reshape2")
library("reshape")
library("lulu")
library("ggrepel")
library("ggh4x")
library("RColorBrewer")
library("rITSx")
library("MicEco")
library("agricolae")
library("ampvis2")
library("pheatmap")
library("tidyverse")

2 Version

1.0

• Support Request Request sent by the user:
  Mona N. Högberg
  

3 Data

96 samples

• Type of data

ITS amplicon

• Data location rackham.uppmax.uu.se /proj/snic2022-22-352

• Uppmax project ID SNIC 2022/22-352

• NGI Project ID P9723

• Database used Unite

4 Tools

NFCore-Ampliseq (Dada2)

#Sample Info
head(sample_info_tab)

5 Workflow

I’m describing here the most important steps in the analysis. Moreover, it is possible to find all parameters and figures at the .rmd file.

5.0.1 Inspect read quality profiles

Quality profiles of the forward reads:

QC Foward Reads.

In gray-scale is a heat map of the frequency of each quality score at each base position. The median quality score at each position is shown by the green line, and the quartiles of the quality score distribution by the orange lines. The red line shows the scaled proportion of reads that extend to at least that position

The forward reads are good quality I truncated the forward reads at position 240.

Quality profile of the reverse reads:

QC Reverse Reads.

The reverse reads are of significantly worse quality, especially at the end, which is common in Illumina sequencing. This isn’t too worrisome, as DADA2 incorporates quality information into its error model which makes the algorithm robust to lower quality sequence, but trimming as the average qualities crash will improve the algorithm’s sensitivity to rare sequence variants. Based on these profiles, I truncated the reverse reads at position 200 where the quality distribution crashes.

5.0.2 Identify primers

The universal primer set 5.8S and ITS1F (Fierer et al., 2005; Yarwood et al., 2010) primers were used to amplify this dataset. We record the DNA sequences, including ambiguous nucleotides, for those primers.

FWD <- “CTTGGTCATTTAGAGGAAGTAA”
REV <- “GCTGCGTTCTTCATCGATGC”

Sanity check of primers and adapters

##                  Forward Complement Reverse RevComp
## FWD.ForwardReads       0          0       0       0
## FWD.ReverseReads       0          0       0       0
## REV.ForwardReads       0          0       0       0
## REV.ReverseReads       0          0       0       0

Success! Primers are no longer detected in the cutadapted reads.

The primer-free sequence files are now ready to be analyzed through the DADA2 pipeline.

5.0.3 Learn the Error Rates

The DADA2 method relies on a parameterized model of substitution errors to distinguish sequencing errors from real biological variation. Because error rates can (and often do) vary substantially between sequencing runs and PCR protocols, the model parameters can be discovered from the data itself using a form of unsupervised learning in which sample inference is alternated with parameter estimation until both are jointly consistent.

In order to verify that the error rates have been reasonably well-estimated, we inspect the fit between the observed error rates (black points) and the fitted error rates (black lines) in the next Figure:

Error Rate Foward Reads.

Error Rate Foward Reads.

Everything looks reasonable and we proceed with confidence.

After filtering, we moved to the next step:

Sequence inference: The DADA2 sequence inference step removed - nearly - all substitution and indel errors from the data. We now merge together the inferred forward and reverse sequences, removing paired sequences that do not perfectly overlap as a final control against residual errors.

5.0.4 Construct sequence table and remove chimeras

The DADA2 method produces a sequence table valued by the number of times each sequence was observed in each sample.

Notably, chimeras have not yet been removed. The error model in the sequence inference algorithm does not include a chimera component, and therefore we expect this sequence table to include many chimeric sequences. We now remove chimeric sequences by comparing each inferred sequence to the others in the table, and removing those that can be reproduced by stitching together two more abundant sequences.

The frequency of chimeric sequences varies substantially from dataset to dataset, and depends on on factors including experimental procedures and sample complexity.

sum(seqtab.nochim)/sum(seqtab)

## [1] 0.950831

5.0.5 Inspect distribution of sequence lengths:

As expected, quite a bit of length variability in the the amplified ITS region.

Distribution of sequence lengths.

5.0.6 Track reads through the pipeline

We now inspect the the number of reads that made it through each step in the pipeline to verify everything worked as expected.

Track Reads.

Looks good! We kept the majority of our raw reads. The most impacting was the filter step.

5.0.7 Taxa prevalence versus total counts.

Taxa prevalence.

Each point is a different taxa. Exploration of the data in this way is often useful for selecting filtering parameters, like the minimum prevalence criteria we will used to filter the data above.

In this study, I removed taxa not seen more than 2 times in at least 1% of the samples. This protects against an ASV with small mean & trivially large C.V.

Count prevalence
- Kingdom = 1
– Phylum = 15
— Class = 53
— Order = 128
—- Family = 232
—– Genus = 458
—— ASVs = 551

5.1 Replicates Evaluation

In this step an evaluation of the replicates has been performed. The main goal was to understand if they are any significant difference between the replicates A and B.

Alpha Diversity Plot comparing the Replicates A and B.

Quality control analysis using matched samples from 3 different Transects: A, B and C of the experiment and replicates samples on each Ecotype. Comparison of alpha diversity in technical replicates samples on all Ecotypes from each Transect. ASV richness and ASV Shannon diversity.

Alpha Diversity.

Evaluation of alpha diversity: looking for difference at the count read distribution

Beta Diversity Plot comparing the Replicates A and B.

Quality control analysis using matched samples from replicates A and B. Beta diversity using Jensen-Shannon distance

Beta Diversity Distance

Beta Diversity Plot comparing the Replicates A and B.

Beta Diversity Replicates.

Beta Diversity Replicates by Type.

In summary: there is no difference between A and B samples. In this way we choose working the replicate A.

Rarefaction

Rarefaction.

Rarefaction.

5.2 Main Analysis

• Alpha Diversity Plot - A.

Raw alpha-diversity, i.e. without rarefying

• comparing alpha diversity based on raw data has the huge problem of ignoring sample size differences. Usually there is a trend that more reads correlate with higher diversity.

• Here, therefore I checked, whether there are significant sample size differences between the groups.

• We further look at the residuals of a linear fit alpha diversity vs sample size to correct for this confounder.

• Check whether sample_sums/sample sizes/library sizes differ between groups

• sample size adjustment has no influence on richness, so sample_sizes are compared for raw counts

AlphaDiversity Main Analysis.

Alpha Diversity = BoxPlot

AlphaDiversity Main Analysis.

summary(aov(df2$Observed ~ Transect * Ecotype, data = df2))

##                  Df Sum Sq Mean Sq F value        Pr(>F)    
## Transect          2   1715     858   0.733      0.488718    
## Ecotype           4 151580   37895  32.400 0.00000000017 ***
## Transect:Ecotype  8  51828    6478   5.539      0.000243 ***
## Residuals        30  35087    1170                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot by mean

AlphaDiversity Main Analysis.

AlphaDiversity Main Analysis.

• Beta Diversity Plot - A.

Beta Diversity Main Analysis.

Beta Diversity Main Analysis.

• Beta Diversity Plot - B.

Beta Diversity Main Analysis.

Ecotype

Beta Diversity Main Analysis.

• Permutational analysis of variance

HeatMap

AlphaDiversity Main Analysis.

AlphaDiversity Main Analysis.

• Relative abundance of the top 20 ASVs

AlphaDiversity Main Analysis.

6 Summary

Help is needed with running the “nfcore/ampliseq” pipeline developed at NGI for the analyses of fungal ITS1 amplicons, Illumina Miseq analysis NGI project ID P5953 (M.Hogberg_17_01_project summary from 2018-01-19 by Chuan Wang refers to P9723, >=Q30 (mean(SD), 70(2) (%), Sum Reads=15 650 000, Mean reads per sample 171 711, 1 pool of amplicons, 1 flowcell v3, PE 2x301bp (validated method), demultiplexing, quality control and raw data delivery on Uppmax (validated method). Agreement number M.Hogberg_16_01-20160826. Grus delivery project delivery 00654. Because my support application 2019-01-17 was rejected, I have collaborated with partners in US on this matter but all is extremely delayed for known reasons. I got some hope today when I read about the recently developed pipeline for fungal analyses! Unfortunately, I have no programming skills but have a BSc in Molecular Biology.

Short summary of the work.

7 Further Work

Further steps to be taken (if needed).

8 References

Relevant references for methods, tools etc.

9 Deliverables

Files delivered to the user with descriptions.

9.1 Directory

/data/processed/b/

8 directories, 18 files

Total size is XX GB.

10 Timeline

11 Practical Info

11.1 Data responsibility

The responsibility for data archiving lies with the PI of the project. We do not offer long-term storage or retrieval of data.

• NBIS & Uppnex: We kindly ask that you remove the files from UPPMAX/UPPNEX. The main storage at UPPNEX is optimized for high-speed and parallel access, which makes it expensive and not the right place for longer time archiving. Please consider others by not taking up the expensive space. Please note that UPPMAX is a resource separate from the Bioinformatics Platform, administered by the Swedish National Infrastructure for Computing (SNIC) and SNIC-specifc project rules apply to all projects hosted at UPPMAX.
  
• Sensitive data : Please note that special considerations may apply to the human-derived legally considered sensitive personal data. These should be handled according to specific laws and regulations as outlined e.g. here.
  
• Long-term backup : We recommend asking your local IT for support with long-term data archiving. Also a newly established Data Office at SciLifeLab may be of help to discuss other options.

11.2 Acknowledgments

If you are presenting the results in a paper, at a workshop or conference, we kindly ask you to acknowledge us.

• NBIS staff are encouraged to be co-authors when this is merited in accordance to the ethical recommendations for authorship, e.g. ICMJE recommendations. If applicable, please include Juliana, Assis Geraldo, National Bioinformatics Infrastructure Sweden, Science for Life Laboratory, NBIS, as co-author. In other cases, NBIS would be grateful if support by us is acknowledged in publications according to this example:

“Support by NBIS (National Bioinformatics Infrastructure Sweden) is gratefully acknowledged.”

• UPPMAX kindly asks you to acknowledge UPPMAX and SNIC. If applicable, please add:

“The computations were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under Project SNIC 2022-22-352.”

• NGI : For publications based on data from NGI Sweden, NGI, SciLifeLab and UPPMAX should be acknowledged like so:

“The authors would like to acknowledge support from Science for Life Laboratory (SciLifeLab), the National Genomics Infrastructure (NGI), and Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) for providing assistance in massive parallel sequencing and computational infrastructure.”

12 Support Completion

You should soon be contacted by one of our managers with a request to close down the project in our internal system and for invoicing matters. If we do not hear from you within 30 days the project will be automatically closed and invoice sent. Again, we would like to remind you about data responsibility and acknowledgements, see sections: Data Responsibility and Acknowledgments.

You are welcome to come back to us with further data analysis request at any time via http://nbis.se/support/support.html.

Thank you for using NBIS.